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Abstract. Consider the problem of optimally matching two measures on the circle, or equiva- 
lently two periodic measures on M, and suppose the cost c(x,y) of matching two points x, y satisfies 
the Monge condition: c(xi, yi) + c(a;2, j/2) < c{xi,y2) + c{x2,yi) whenever xi < X2 and yi < y2. We 
introduce a notion of locally optimal transport plan, motivated by the weak KAM (Aubry— Mather) 
theory, and show that all locally optimal transport plans are conjugate to shifts and that the cost of 
a locally optimal transport plan is a convex function of a shift parameter. 

This theory is applied to a transportation problem arising in image processing: for two sets of 
point masses on the circle, both of which have the same total mass, find an optimal transport plan 
with respect to a given cost function c satisfying the Monge condition. In the circular case the sorting 
strategy fails to provide a unique candidate solution and a naive approach requires a quadratic number 
of operations. For the case of N real-valued point masses we present an 0{N\ loge|) algorithm that 
approximates the optimal cost within e; when all masses are integer multiples of 1/Af, the algorithm 
gives an exact solution in 0{N log M) operations. 
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1. Introduction. The transport optimization problem, introduced by G. Monge 
in 1781 and shown by L. Kantorovich in 1942 to be an instance of linear programming, 
is a convex optimization problem with strong geometric features. A typical example is 
minimization of mean-square displacement between two given finite marginal measures 
supported on convex compacta in Euclidean space: in this case a solution is defined 
by gradient of a convex function that satisfies a suitable Monge-Ampere equation. 
Various generalizations of this result and rich bibliographies can be found, e.g., in the 
survey |12j or the recent monograph 19 . 

Further constraints on the two marginals or their supports may furnish the prob- 
lem with useful additional convex structure. One way of making this statement 
quantitative is to consider the algorithmic complexity of the corresponding numer- 
ical transport optimization schemes. In particular when the two measures live on 
segments of straight lines, the optimal map is monotone and may be found by sort- 
ing, which takes O(nlogn) operations when the data come in the form of discrete 
n-point histograms. If the input data are already sorted, this count falls to 0{n). 

The optimal transport is well understood also when the marginals live on a com- 
pact Riemannian manifold [lOj ; the existence and characterization of optimal map in 
the case of a flat torus and quadratic cost have been established a decade ago [8]. 
However, the algorithmics of even the simplest setting of the unit circle is no longer 
trivial, because the support of the measures is now oriented rather than ordered. A 
naive approach would require solving the problem for each of n different alignments 
of two n-point histograms, thereby involving OirP') operations. 
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In this paper we present an efficient algorithm of transport optimization on the 
circle, which is based on a novel analogy with the Aubry-Mather (weak KAM) theory 
in Lagrangian dynamics (see, e.g., [3J HI HI])- The key step is to lift the transport 
problem to the universal cover of the unit circle, rendering the marginals periodic and 
the cost of transport infinite. However, it still makes sense to look for those transport 
maps whose cost cannot be decreased by any local modification. Different locally 
optimal maps, which cannot be deformed into each other by any local rearrange- 
ment, form a family parameterized with an analogue of the rotation number in the 
Aubry-Mather theory. One can introduce a counterpart of the average Lagrangian, 
or Mather's a function in the Aubry-Mather theory, which turns out to be efficiently 
computable. As we show below, its minimization provides an efficient algorithm of 
transport optimization on the circle. The class of cost functions for which this theory 
works includes all costs with the Monge property, such as the quadratic cost or costs 
generated by natural Lagrangians with time-periodic potentials [5J [Tl] . 

Note that the problem of optimally matching circular distributions appears in a 
variety of applications. Important examples are provided by image processing and 
computer vision: image matching techniques for retrieval, classification, or stitching 
purposes [2TJ |6] are often based on matching or clustering "descriptors" of local fea- 
tures [15) , which typically consist of one or multiple histograms of gradient orientation. 
Similar issues arise in object pose estimation and pattern recognition [15l[Tl]. Circular 
distributions also appear in a quite different context of analysis of color images, where 
hue is parameterized by polar angle. In all these applications, matching techniques 
must be robust to data quantization and noise and computationally effective, which 
is especially important with modern large image collections. These requirements are 
satisfied by the optimal value of a transport cost for a suitable cost function. 

This paper is organized as follows. In ij2] we give a specific but nontechnical 
overview of our results. After the basic definitions are given in Sj3l including that of 
locally optimal transport plans, in f|4]we give an explicit description of the family of 
locally optimal transport plans: they are conjugate, in measure theoretic sense, to 
rotations of the unit circle (or equivalently to shifts of its universal cover) . This result 
is in direct analogy with conjugacy to rotations in the one-dimensional Aubry-Mather 
theory |3] . As shown in ^ the average cost of a locally optimal transport plan is a 
convex function of the shift parameter. Moreover, the values of this function and its 
derivative are efficiently computable when the marginal measures are discrete, which 
enables us to present in [jHl a fast algorithm for transport optimization on the circle. 
The same section contains results of a few numerical experiments. Finally a review 
of related work in the computer science literature is given in SjT] 

2. Informal overview. For two probability measures /io, Ai on the unit circle 
T = R/Z and a given cost c{x,y) of transporting a unit mass from i to y in T, the 
transport cost is defined as the inf of the quantity 



over the set of all couplings 7 of the probability measures /to, /ti (i.e., all measures on 
TxT with marginals Ao, Ai)- These couplings are usually called transport plans. 

Suppose that the cost function c(-, •) on T x T is determined via the relation 
c{x, y) — inf c(x, y) by a function c(-, •) on RxR satisfying the condition c{x+l, y+1) = 
c{x, y) for all x, y; here inf is taken over all x, y whose projections to the unit circle 
coincide with x, y. We lift the measures Ao and Ai to M, obtaining periodic locally 




(2.1) 
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Fig. 2.1. Construction of the locally optimal transport plan 79, 

finite measures /io, /ii, and redefine 7 to be their coupling on M x M. It is then 
convenient to replace the problem of minimizing the integral (|2.ip with "minimization" 
of an integral 

^(7)= // c{x,y)j{dxxdy). (2.2) 
J Jrxs. 

Although the latter integral is infinite, it still makes sense to look for transport plans 7 
minimizing / with respect to local modifications, i.e., to require that for any compactly 
supported signed measure S of zero mass and finite total variation, the difference 
7(7 + S) — 7(7) , which is defined by a finite integral, be nonnegative. These locally 
optimal transport plans are the main object of this paper. 

Assume that the cost function c(x, y) satisfies the Monge condition (alternatively 
known as the continuous Mongc property, see [Hl^)" 

c(xi,yi) + c{x2,y2) < c{xi,y2) + c{x2,yi) (2.3) 

for all xi < X2 and yi < y2- An example of such a cost function is \x — y\^, where 
A > 1; in this case the quantity MKx{fiQ, fti) — (inf^ 1(7))^/'*' turns out to be a metric 
on the set of measures on the circle, referred to as the Monge-Kantorovich distance 
of order A. The value A = 1 can still be treated in the same framework as the limiting 
case A 1; it is sometimes called the Kantorovich-Rubinshtein metric or, in image 
processing literature, the Earth Mover's distance [TB]. 

The Monge condition (12. 3p implies that whenever under a transport plan the 
mutual order of any two elements of mass is reversed, the transport cost can be 
strictly reduced by exchanging their destinations. It follows that a locally minimal 
transport plan moves elements of mass monotonically, preserving their spatial order. 

The whole set of locally optimal transport plans for a given pair of marginals 
/io, /ii can be conveniently described using a construction represented in fig. 12.11 
Let Fq, Fi be cumulative distribution functions of the measures /io, /ii normalized 
so that Fo{0) ~ Fi{0) ~ 0. We shall regard graphs of Fq, Fi as continuous curves 
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including, where necessary, the vertical segments corresponding to jumps of these 
functions (which are caused by atoms of /iq, fJ-i)- Each of these curves specifies a 
correspondence, Fq^ or F-^^, between points of the vertical axis Ov, representing 
elements of mass, and points of the horizontal axis Ou, representing spatial locations, 
and maps the Lebesgue measure on Ov into fiQ or fii on Ou. This correspondence is 
monotone and defined everywhere except on an (at most countable) set of v values 
that correspond to vacua of the measure in the Ou axis. 

Define now F^{u) = Fi{u) — 9. Then {Ff)~^ represents a shift of the Ov axis by 9 
followed by an application of the correspondence F^^, and still induces on the Ou 
axis the same measure fii as Fi. A transport plan that takes an element of mass 
represented by v from F^^{v) to {Ff )~^{v) is, by construction, a monotone coupling 
of fj,Q and /ii, and thus a locally optimal transport plan. Moreover, it is shown in ^ 
that all locally optimal transport plans can be obtained using this construction as the 
parameter 9 runs over M. 

Finally define the average cost C[Pg^p^]{9) of the plan per unit period: 

ClFo,F,m = f' c{F,-\v),{Flr\v))dv. 
Jo 

It is shown in Sj5]that the Mongc condition implies convexity of C[pg p-^]{9) and that 
the global minimum of this function in 9 coincides with the minimum value of the 
transport cost on the unit circle (|2.ip . 

When the marginals /ip, Mi are purely atomic with finite numbers Uq and ni of 
atoms in each period, the function C becomes piecewise affine. In fj^l we present an 
algorithm to approximate its minimum value to accuracy e, using a binary search that 
takes 0{{nQ + ni)log(l/e)) operations in the real number computing model. When 
masses of all atoms are rational numbers with the least common denominator M, 
this search returns an exact solution provided that e < 1/M. This gives an 0{{no + 
ni) logM) exact transport optimization algorithm on the circle. 

3. Preliminaries. Let T = R/Z be the unit circle, i.e., the segment [0, 1] with 
identified endpoints. By tt: K T denote the projection that takes points of the 
universal cover R to points of T. 

3.1. The cost function. A cost function is a real-valued function c(-, ■) defined 
on the universal cover R of the circle T. We assume that it satisfies the Monge 
condition: for any xi < X2 and j/i < 2/2, 

c{xi,yi) + c{x2,y2) - c{xi,y2) - c{x2,yi) < 0. (3.1) 

Additionally c is assumed to be lower semicontinuous, to be invariant with respect to 
integer shifts, i.e., 

c{x + l,y + l)^c{x,y) (3.2) 

for all X, y, and to grow uniformly as |a; — y| — t- oo; for any P there exists a fi- 
nite R{P) > such that 

c{x,y) > P whenever \x — y\ > R{P)- (3-3) 

Note that the latter condition implies that the lower semicontinuous function c is 
bounded from below (and guarantees that the minima in a number of formulas below 
are attained). 
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Note that the Monge condition (|3.ip holds for any twice continuously differentiable 
function c such that d^c{x, y)/dxdy < 0. If the cost function depends only on x — y, 
this reduces to a convexity condition: —d^c{x — y)/dxdy = c"{x — y) > 0. In 
particular, all the above conditions are satisfied for the function c{x,y) ~ \x — y\^, 
which appears in the definition of the Monge-Kantorovich distance with A > 1, and, 
more generally, for any function of the form c{x — y) + f (x) + g{y) with strictly convex c 
and periodic / and g. 

For a cost function c satisfying all the above conditions, the cost of transporting 
a unit mass from a; to y on the circle is defined as c{x,y) — vaic{x,y), where x, y 
are points of T and inf is taken over all x, y in M such that irx = x and iry = y. 
Using the integer shift invariance, this definition can be lifted to the universal cover 
as c{x, y) = inffcgz c(x, y + k). 

Condition is all that is needed in §31 which is concerned with locally optimal 
transport plans on R. Conditions p.2p . (|3.3|) come into play in fJSJ which deals with 
transport optimization on the circle. 

3.2. Distribution functions. For a given locally finite measure /i on M define 
its distribution function F^j, by 

F^(0) = 0, F^{x)^ ti{{0,x]) ioT x>0, F,,{x)^-ti{{x,0])foT x<0. (3.4) 

Then iJ,{{xi,X2]) — F^(x2) — F^(xi) whenever xi < X2, and this identity also holds for 
any function that differs from F^ by an additive constant (the normalization F^ (0) = 
is arbitrary). When fi is periodic with unit mass in each period, it follows that for all 
a; in M 

F^{x + l)=F^{x) + l. (3.5) 

The inverse of a distribution function i^^ is defined by 

F-\y) = inf{a;: y < F^{x)} = sup{x: y > F^ix)}. (3.6) 

Definitions (13. 4p and p.6p mean that i^^, F^^ are right-continuous. Discontinuities 
of Ff^ correspond to atoms of /i and discontinuities of its inverse, to "vacua" of /i, i.e., 
to intervals of zero fi measure. 

For a distribution function F^ define its complete graph to be the continuous curve 
formed by the union of the graph of Ff^ and the vertical segments corresponding to 
jumps of F^. Accordingly, by a slight abuse of notation let F^{{x}) denote the set 
[F^{x - 0),F^{x)] (warning: F^{{x}) D {F^{x)}) and let F^{A) = [j.eAF^ii^}) ioi 
any set A. 

3.3. Local properties of transport plans. Let fiQ, fii be two finite positive 
measures of unit total mass on T and /ig, /ii their liftings to the universal cover M, 
i.e., periodic measures such that fJ.i{A) = fli{TTA), i = 0, 1, for any Borel set A that 
fits inside one period. Periodicity of measures here means that fi(A + n) — fi{A) for 
any integer n and any Borel A, where A + n = {x + n: x G A}. 

Definition 3.1. A (locally finite jl] transport plan with marginals /io and /ii is 
a locally finite measure 7 on R x R such that 

(i) for any x m R the supports of measures 7((— oo,a:] x •) and "/{■ x (—oo,x]) 
are bounded from above and the supports of measures "f{{x,oo) x •), "/{■ x (a;, 00)) are 
bounded from below; 



In what follows the words 'locally finite' defining a transport plan will often be dropped. 
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(ii) j{A X M) = ^o(^) and j(R x B) ^ fJ-i{B) for any Borel sets A, B. 

The quantity ^{A x B) is the amount of mass transferred from A to i? under 
the transport plan 7. Condition ([l} imphes that the mass supported on any bounded 
interval gets redistributed over a bounded set (indeed, a bounded interval is the 
intersection of two half-lines), but is somewhat stronger. 

Definition 3.2. A local modification 0/ the locally finite transport plan j is a 
transport plan 7' such that 7 and 7' have the same marginals and j' — j is a compactly 
supported finite signed measure. A local modification is called cost-reducing if 

J J c{x, y) (7'(da:: x dy) - 7(da:: x dy)) < 0. 

A locally finite transport plan 7 is said to be locally optimal with respect to the cost 
function c or c-locally optimal if it has no cost-reducing local modifications. 

4. Conjugate transport plans and shifts. Let Ua, Ui be two copies of R 
equipped with positive periodic measures /io, Mi whose distribution functions Fq, Fi 
satisfy p.Sp . so that all intervals of unit length have unit mass. Let furthermore Vb, Vi 
be two other copies of R equipped with the uniform (Lebesgue) measure. 

4.1. Normal plans and conjugation. We introduce the following terminology: 
Definition 4.1. A locally finite transport plan v on Vq x Vi with uniform 

marginals is called normal. 

Definition 4.2. For a normal transport plan v its conjugate transport plan 

y{Fo,Fx\ jg ^ transport plan on Uq x Ui such that for any Borel sets A, B 

„lF'o.F^]^ X B) !^{Fo{A) X Fi{B)). (4.1) 

Lemma 4.3. For a normal transport plan v its conjugate IS a locally finite 

transport plan on Uq x Ui with marginals ^0, ^i. 

Proof. Since distribution functions _Fo, Fi and their inverses preserve bound- 
edness, condition ([l} of Definition 13.11 is fulfilled. Definition 14.21 condition ^ of 
Definition 13.11 and formula (|3.4p together imply that 

i^[^°'^il((«i,M2] X Ui) = iy{Fo{{ui,U2]) X Fi([/i)) = i^{[Fo{ui), Foiu^)] x ^i) 
= -^0(^2) - Fa{ui) = ^J.a{{ul,U2]). 

Similarly v^^O'^^^^Uq x (mi,W2]) = fJ.i{{ui,U2]). Thus i^i^"^'^^^ satisfies condition ^ of 
Definition 13.11 on intervals and therefore on all Borel sets. □ 

Lemma 4.4. For any transport plan j onUo x Ui with marginals /io o,n,d /ii there 
exists a normal transport plan v such that 7 is conjugate to v: 7 = i/^^o-'^^^. 

Proof. For non-atomic measures /xq and fj,i the required transport plan is given 
by the formula iy{A x B) ^ ^{Fq^{A) x F{~^{B)), which is dual to (|44|) . However if, 
e.g., /io has an atom, then the function Fq^ is constant over a certain interval and 
maps any subset A of this interval into one point of fixed positive measure in Uq , so 
information on the true Lebesgue measure of A is lost. In this case extra care has to 
be taken. 

Recall that a locally finite measure has at most a countable set of atoms. Let 
atoms of /io be located in (0, 1] at points ui, U2, . ■ . with masses mi, m2, .... Since 
7({ui} X Ui) = ^o{{ui}) = TUi > 0, there exists a conditional probability measure 
pi' I — 7({'"i} ^ ■)/n^i- For a set A C (0, 1] define a "residue" transport plan 

7(A X S) - 7(A xB)-J2^ m, 5^, {A) p{B \ u,), 
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where (5„ is the Dirac unit mass measure on Uq concentrated at u, and extend 7 to 
general A using periodicity. We thus remove from 7 the part of 7 whose projection 
to the first factor is atomic. Define a transport plan k on Vq x ?7i by 

k{C xB)^ HC n Fo({u,})) p{B I u,) + j{F^-\C) x B), 

where C is a Borel set in Vq and A(-) denotes the Lebesgue measure in Vq. Clearly 
k{Fo{A) xB) = 7(A X B). Repeating this construction for the second factor, with k in 
place of 7, we get a normal transport plan v such that j{A x B) — i'{Fq[A) x Fi[B)). 
□ 

Since we are ultimately interested in transport optimization with marginals /xq , Mi 
rather than with uniform marginals, two normal transport plans vi, 1/2 will be called 
equivalent if they have the same conjugate. Two different normal transport plans 
can only be equivalent if one or both measures /ig or /.ti have atoms, causing loss 
of information on the structure of v in segments corresponding to these atoms. The 
proof of Lemma 14.41 gives a specific representative of this equivalence class of normal 
plans. 

4.2. Locally optimal normal transport plans are shifts. Fix a cost function 
c: J7o X [/i — > R that satisfies the Monge condition p.ip and define 

c^FoM^o.v,) = c{F,-\vo),F,-'{v,)). (4.2) 

For non-atomic measures fio, fii, it satisfies the Monge condition 

C[Fo,F^]iv' ,w') + C[p^,^p^]iv" ,w") - C[f^,^f^]{v' , w") - C[Fo,F^]{v" < 

whenever v' < v" and w' < w"; this inequality can only turn into equality if either 
v' , v" or w' , w" correspond to an atom of the respective marginal (/iq or /ii) of i/I^f'^i] , 
i.e., if C[Fo,Fi] is constant in either first or second argument. In spite of this slight 
violation of definition of §3.11 we will still call cyFo,Fi\ a cost function. 

Here and below, variables w, u' , uq, ui, . . . are assumed to take values in Uq or Ui 
and variables v, v' , vq, vi, . . . , w, w' , . . . , in Vq or Pi. 

Lemma 4.5. A transport plan ^ on Uq x Ui with marginals fiQ, fii is c-locally 
optimal if and only if it is conjugate to a C[FQ.Fi]-locally optimal normal transport plan 
V. In particular, all normal transport plans with the same locally optimal conjugate 
are locally optimal. 

Proof. Note that v' — v \s compactly supported if and only if the difference of 
the respective conjugates 7' — 7 is compactly supported. The rest of the proof follows 
from the identity 

j j C[Fq,Fi]{viiV2) {v'i'ivi X dv2) ~ viAvi x dv2)) 

= j j c(ui, U2) (7'(dui X du2) - 7(dMi x du2)) 

established by the change of variables vi — Fq{ui), V2 = -^1(^2) (here jumps of the 
distribution functions are harmless because C[Fo,Fi] is constant over respective ranges 
of its variables) . □ 

Transport optimization with marginals ^o, fJ-i is thus reduced to a conjugate 
problem involving uniform marginals and the cost C[Fo,Fi]- It turns out that any 
C[Fo,Fi]-0'ptinial normal transport plan must be supported on a graph of a monotone 
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function, and due to uniformity of marginals this function can only be a shift by a 
suitable real increment 0. More precisely, the following holds: 

Theorem 4.6. Let no, jii he two periodic positive measures defined respectively 
on Uq, Ui with unit mass in each period and let Fi: Ui — > Vi, i = 0,1, be their 
distribution functions. Then any C[j^q j?^] -locally optimal normal transport plan on Vq x 
Vi is equivalent to a normal transport plan vg with supp vq = {{v,w): w = v + O}, and 
conversely vg is c^p^ p^-^-locally optimal for any real 9. All c-locally optimal transport 
plans on Uq x Ui with marginals /ig, Hi are of the form — {vg)^-^"'-^'-^ . 

The proof, divided into a series of lemmas, is based on the classical argument: a 
nonoptimal transport plan can be modified by "swapping" pieces of mass to render its 
support monotone while decreasing its cost. This argument, carried out for plans with 
uniform marginals on Vq x Vi , is combined with the observation that a monotonicaly 
supported plan with uniform marginals can only be a shift. Then Lemma 14.51 is used 
to extend this result to transport plans on Uq x Ui. 

Throughout the proof fix a normal transport plan v and define on Vq x Vi the 
functions 

ri,(v,w) = oo,w] X (w, oo)), l^{v,w) = h'{{v,oo) x (— oo,it;]). (4-3) 

To explain the notation r,^, /„m observe that, e.g., rj^(u, w) is the amount of mass that 
is located initially to the left of v and goes to the right of w. 

Lemma 4.7. The function r^, (resp. lu) is continuous and monotonically increas- 
ing in its first (second) argument and is continuous and monotonically decreasing in 
its second (first) argument, while the other argument is kept fixed. 

Proof. Monotonicity is obvious from (j4.3|) . To prove continuity observe that 
the second marginal of v is uniform, which together with positivity of all involved 
measures implies that in the decomposition 

v{Vo X • ) = iy{{-oo, w] X • ) + iy{{v, oo) x • ), 

both measures in the right-hand side cannot have atoms. This implies continuity 
of r^, li, with respect to the second argument. A similar proof holds for the first 
argument. □ 

Lemma 4.8. For any v there exist Wu{v) and mu{v) > .such that 

r„{v, w„{v)) = li,{v, Wi,{v)) = mi,{v). (4.4) 

The correspondence v t-> w,^{v) is monotone: Wiy{vi) < Wi,{v2) for vi < V2. 

Proof. Clearly r^{v, —oo) ~ oo, r-^(w, oo) — 0, li,{v, — oo) — 0, l^{v, oo) — oo. The 
continuity of the functions r^{v, •), /^(w, •) in the second argument for a fixed v implies 
that their graphs intersect at some point {wy{v),mi,{v)), which satisfies (|4.4p . Should 
the equality r^{v,w) — li,{v,w) hold on a segment [w',w"], we set w,y{v) to its left 
endpoint w'; this situation, however, will be ruled out by the corollary to Lemma 14.101 
below. Monotonicity of Wu{v) follows from monotonicity of ru{-,w), l^{-,'w) in the 
first argument for a fixed w. indeed, for V2 > vi the equality r^(v2,w) — l,y{v2,w) is 
impossible for w < Wu{vi) because for such w we have r^{v2,w) > r^{vi,w^{v)) = 

liyivi,W,,{v)) >li,{v2,w). □ 

Equalities (j4.4p mean that the same amount of mass m^{v) goes under the plan v 
from the left of v to the right of Wv{v) and from the right of v to the left of w^iv). 
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We are now in position to use the Monge condition and show that this amount can 
be reduced to zero by modifying the transport plan locally without a cost increase. 

Lemma 4.9. For any v there exists a local modification Vy of v such that w^^ (v) = 
w^{v) (with Wi, defined as in Lemma \4-.8\ l, m^^{v) = 0, and Vy is either cost-reducing 
in the sense of Definition \3.2\ or is equivalent to v. 

Proof. Let w — w^{v) and m = m^{v). If m = 0, there is nothing to prove. 
Suppose that m > and define 

— sup{w' : lu{v, w') — 0}, u>+ = inf{?«' : r^{v, w') — 0}, 
v~ = sup{u' : r^{v' , w) = 0}, = inf{u' : l^{v' , w) = 0}. 

By local finiteness of the transport plan v all these quantities arc finite. Since to > 0, 
continuity of r^, 1,^ implies that the inequalities < w < and < v < v'^ are 
strict. Consider the measures 

{■) = h'{ ■ X (w, w^)) on P'^i') = x {w~ ,w)) on 

cr^ (•) — iy{{v,v'^) X ■ ) on w), ~ iy((v^ ,v) x • ) on {w, w^). 

Equalities (14. 4p mean that all these measures have the same positive total mass m. 
Note that the Lebesgue measures of intervals {v,v~^), {w~,w), and {w,w'^) 

may be greater than to, because some mass in these intervals may come from or go 
to elsewhere. 

The functions ?"„,(•) ~ ri,{-,w), ly{-) = luiv,') are monotonically increasing and 
fv{ ) = i'i,{v,-), lu]{-) ^ li,{-,w) are monotonically decreasing, with their inverses r~^, 
l^^ , r^^, l~^ defined everywhere except on an at most countable set of points. These 
functions may be regarded as a kind of distribution functions for the measures p~ , 
fj^, (T+, p"*" respectively, mapping them to the Lebesgue measure on (0,to). 

Under the plan v, mass to is sent from (v~,v) to {w,w~^) and from (v,v~^) to 
{w~,w). We now construct a local modification Vy of the transport plan v that 
moves mass to from the interval {v^,v) to {w^,w) and from to {w,w^), and 

show that it is cost-reducing unless measures po , pi have atoms corresponding to the 
intervals under consideration. 

Observe first that the normal plan ly induces two transport plans , t; that map 
measures p~ to cr+ and p+ to correspondingly: 

Tr{A X B) = v{A n (-oo,t;) x Br\{w, +oo)) = v{A n {v~ ,v) x Bn (w, 
Ti{A X B) = v{Ar\ (i;,+oo) x Bn (-oo,w)) = v{Ar\ {v,v+) x B f] 

where A C B C (w^,^^) are two arbitrary Borel sets and the n opera- 

tion takes precedence over x . By an argument similar to the proof of Lemma 14.41 
there exist two transport plans Xr and xi mapping the Lebesgue measure on (0, to) 
respectively to cr^, and such that 

TriA X B)= Xr{rwiAn{v~,v)) X Bn{w,w+)), 
Ti{A xB)^ xi{lw{Ar\(v,v+)) xBn{w-,w)). 

Define now two transport plans f/ , fv that send mass elements to the same destinations 
but from interchanged origins: 

fr{A xB)= Xr{lw{A n {v,v+)) X B Ci {w , w+)), 
fi{A X B) = xi{rvj{Ar\{v-,v)) x Bn{w-,w)). 
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This enables us to define 

I'yiA X B)^ v{A X B)- Tr{A X B) - Ti{A X B) + fr{A x B)+fi{Ax B). 

Since rr(AxR) = fr{AxW) = p~{A) etc., the transport plan has the same uniform 
marginals as v, i.e., it is a local modification of v. Observe furthermore that by the 
construction of no mass is moved under this plan from the left-hand side of v to 
the right-hand side of w and inversely, i.e., that m^^(v) = 0. 

It remains to show that v^, is either a cost-reducing modification of v or equivalent 
to it. By the disintegration lemma (see, e.g., [2]) we can write Xr{da x dw') = 
da dGr{w' I a) and x/(da x dw') = da dG';(w' | a), where Gr{- \ a) (resp. G;( • | a)) 
are distribution functions of probability measures defined on [w,w"'"] (resp. [w~,w]) 
for almost all < a < m. Denote their respective inverses by G~^{- \ a), G'f^{- \ a) 
and observe that w- < G^\l3' | a) < w < G-i(/3" | a) < w+ for any /?', (3". Thus 



c{v',w') Tr{dv' X dw') ~ 1 1 c{r^^{a),w') Xr(da x dw') 

da / c{r-\a),w')dGr{w' \ a) 

da C dp c(r-i(a),G;i(/3|a)), 
Jo 

where we write c instead of C[p^^p-^\ to lighten notation, and similarly 

JJ c{v',w')Ti{dv' xdw')^ da dp c{C{a),Gi\l3\a)), 

c{v',w') fr{dv' xdw')= da d/S c{l-\a),G;\l3 \ a)), 

Jo Jq 

c{v',w')fi{dv' xdw')^ da dp c{r-\a),Gj-^{p\a)). 
Jo Jo 

The integral in Definition 13.21 now takes the form 

c{v' , w') (lyyidv' X dw') - i^idv' x dw')) 

c(u',w') (^—Tr{dv' X dw') — Ti{dv' X dw') 
+ fr{dv' X dw') -I- fi{dv' X dw')) 

n pi 

da / dp {~c{r-\a),G~\p \ a)) - c{lz\a),G;\p \ a)) 
Jo 

+ c{l-\a),G;\p I a)) + c{r-^{a),GY\p \ a))). 

As r-i(a) <v< l~^{a) and Gf ^(,3 | a) < w < G~^{P \ a) for all a, P, the Monge 
condition (j3.ip implies that cither the value of this integral is negative or the function 
c (i.e., C[pg p-^]) is constant in at least one of its arguments. In the former case the 
transport plan Vy is a cost-reducing local modification of v. in the latter case Vy is 
equivalent to j/. □ 

Lemma 4.10. For any v' < v" there exists a local modification v^'y of v such 
that Wn I ,„ (w) = w^(v) for v' < v < v" , , ,,(v') — rUy , ,,{v") — 0, and in 
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the strip v' < v < v" the support of Vyi coincides with the complete graph of the 
monotone function u'^(-). 

Proof. Let {vi} be a dense countable subset of [v',v"] including its endpoints. 
Set vq = v and define Vi recursively to be the local modification of Vi-i given by 
the previous lemma and such that Wi,.{vi) — Wi,[vi) and m,y.{vi) — 0. Then all Vi 
are either cost-reducing or equivalent to ly and w,y-(vi) — w,y{vi), m^.{vi) — for 
all j > i. Indeed, denote Wi — Wi^-^Vi) and observe that if e.g. Vj > Vi, then, as 
'mi,-{v.i) — r^.{vi,'Wi) — 0, mass from (—00,1;^] does not appear to the right of Wi 
and so does not contribute to the balance of mass around Wj. Therefore for any 
j the possible modification of I'j-i is local to the interval {vi'^Vi"), where Vi' = 
max{ui : i < j,Vi < Vj} and Ui" = m\Ti{vi: i < j,Vi > Vj} (with max and min of 
empty set defined as v' and v"). Thus there is a well-defined limit normal transport 
plan Vac that is either a cost-reducing local modification or equivalent to v and is 
such that, by continuity of the functions r^^ and l^^ in the first argument, (v) 
vanishes everywhere on [v\v"]. 

Consider now the function w^^{-), which coincides with w^{-) on a dense subset 
of so that their complete graphs coincide. For any quadrant of the form 

(— oo,wo) X (u;o,oo) such that wq > Wi,^(uo)i monotonicity of r^^ in the second 
argument implies that 

< Voo{{-co,vo) X {wo,co)) = r^^{vo,wo) < r^^{vo,w^^{vo)) = 0, 

i.e., z/oo((— 00, uq) X (100,00)) = 0. Similarly z^oo((wo,oo) x (—00,^0)) = for any 
quadrant with wq < w^^{vo). The union of all such quadrants is the complement 
of the complete graph of the function v 1— >■ w^{v); this implies that Voo is supported 
thereon. □ 

Corollary 4.11. For any normal transport plan v there exists a real number 
Oi, such that w^{v) = v + 9^. 

Proof. It is enough to show that Wi,{v') ~ v' — w^{v") — v" for all tj', v" . Let 
v' < v" and I'y'^y" be the local modification constructed in the previous lemma. 
Since it has uniform marginals and monotone support, we have w^{v") — w^{v') = 
t^v',v"{{v\v") X {w^{v'),w„{v"))) = v" — v', which completes the proof. □ 

We call the parameter 9i, the rotation number of the normal transport plan v. 

Definition 4.12. A normal transport plan consisting of a uniform measure 
supported on the line {{v,w): w — v + 0} is called a shift and denoted by vg. 

Lemma 4.13. For any 9 the shift vg is c^p^ p-^ylocally optimal. 

Proof. Let v he a, local modification of vg such that the signed measure vg — v is 
supported in (w', v") x (w', w"). Let Vy' y be a local modification of v constructed in 
Lemma 14.101 it coincides with ug over v' < v < v" , and hence everywhere. Since it is 
either cost-reducing or equivalent to 9, it follows that D cannot be cost-reducing with 
respect to vg, i.e., that vg is a cost minimizer with respect to local modifications. □ 

Lemma 4.14. Any C[p^^,p-^ylocally optimal normal transport v with rotation num- 
ber 9 = 9u is equivalent to the shift vg. 

Proof. Let — —v[ = i for j = 1, 2, . . . . All local modifications Vi — Vy'y of v 
constructed as in Lemma 14.101 cannot be cost-reducing and are therefore equivalent 
to v. On the other hand, this sequence stabilizes to the shift vg on any bounded 
subset of Vq X Vi as soon as this set is covered by {—i,i) x {—i,i). Therefore vg has 
the same conjugate as all Vi and is equivalent to i^. □ 

Lemmas 14.131 14.141 and 14.51 together imply Theorem 14.61 



12 



J. DELON, J. SALOMON AND A. SOBOLEVSKI 



5. Transport optimization for periodic measures. Let now c be a cost 
function that satisfies the Monge condition p.ip . the integer shift invariance condi- 
tion p.2p . the growth condition l\'3.'6\i . and is bounded from below. Suppose that 7^ 
is a locally optimal transport plan on Uq x Ui with marginals /xq, fJ-i conjugate to the 
shift h'e- Define C[p^^p-^] as in (|4.2p and let Ff{u) = Fi{u) — 9 as illustrated in fig. 12.11 

Definition 5.1. We call the quantity 

ClFoM(^) - f c^Fo,F,]{v\v' + e)Av' ^ f\{F,-\v'),{F^)-\v'))dv' (5.1) 



the average cost (per period) of the transport plan 79. Observe that it is indifferent 
whether to integrate here from to 1 or from v to v + 1 for any real v. Examples 
of average cost functions C^p^^p^^ for different marginals /io, Mi and different cost 
functions c are shown in fig. 15.11 

The following technical lemma provides a "bracket" for the global minimum 
of C[pg^p^] and estimates of its derivatives independent of /io, /^i- 

Lemma 5.2. The average cost C[pg^p-^] is a convex function that satisfies the 
inequalities 

inf c(x, y) < Ci9) < C[p„,p,] {9) < C{9) (5.2) 

with 

C_{9) = inf c(mi,W2), C{9) = sup c(mi,W2). (5.3) 

-l<ui<2 -l<ui<2 

There exist constants Q_ < Q and L,L>0 such that the global minimum of C is 
achieved on the interval Q] and 

-L< C[p^,^p,^ (6 - 0) < < (9 + 0) < I, (5.4) 

where C[p^^ p{^{') I'he derivative of C^p^^^p-^^ These constants are independent on 
^0, /ii and are given explicitly by formulas (j5.5p . (|5.6p and (j5.7p below. 

The bounds given in the present lemma are rather loose. E.g., for c(x, y) = |x— yj" 
with a > 1, they are C(6') = convmin(|6' + 3|", |6'-3|"), C(6l) = max(|6l + 3|", |6»-3|"), 
and —Q = Q = 6. For symmetric costs like this one it is often possible to replace 
[9, 9] by the interval [—1,1] which may be tighter. 

Proof. To prove convexity of C[pg^p^] it is sufficient to show that CiPg^p^](^{9' + 

0")) < i(qFo,Fi](f') +C[Fo,Fi](6'")) for all 9', 9". Let 9' < 9\ denote 9 =\{9' + 9") 
and write 

C[Fo,Fi](^') = / C[Fo^Fi](w,W + 6')dw = / C[p„^p^^{v' ,v' + 9)<lv' , 



ClFoM0')= C[p^^p^]{v\v' + 9')dv', qFo,F,](n=/ C[p^,P,]{v,V + 9")dv 

Making the change of variables v' — v + 9 — 9' and taking into account that 9 — 9' + 9 = 
29-9' = 9", we get 

2C[Fo,Fi](^) - C'[Fo,Fi](^'') - C![Fa,Fi]iO") 

(c[Fo,Fi](w, V + 9)+ C[F(,,Fi](« + 9-9',V + 9") 

- C[Fo,Fi](w + 9 - 9' ,V + 9) - C[Fo,Fi](w, V + 9")) dv. 
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■ , , ■ , ■ , , , ■ 1 Q^l , , , , , , , , , 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 

(a) Black bars: fi'g, light bars: fi'^ (b) C[j?/ for c{x,y) = \x — y\ 




■-1 -0.5 0.5 1 1.5 2 

(e) C[^^/ J?//] (•) for a non-symmetric cost 

c{x,y) = |0.5 + x- y\^'^ + 0.1 cos(27rx + 1) - 
0.3 sin(27rj/ - 0.5) 

Fig. 5.1. Average cost functions C^p^ p^] for some atomic marginals and costs. The cost \x — y\ 
is regarded as lim \x — y]-^ , A > 1, as A — >■ 1. 
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Since v + 6 — 9' > v and v + 9" > v + 9, the Monge condition for c implies that the 
integrand here is negative on a set of nonzero measure, yielding the desired inequality 
for the function C[p^^p-^]. Note that convexity of C[Fo,Fi] implies its continuity because 
C[Fo,Fi] is finite everywhere. 

Bounds (|5.2p on C[Pg^p^]{9) follow from (|5.ip with v — because v' — 1 < 
Fo^{v') <v' + l,v' + 9-l< {Ff)-^{v') < w' + 6* + 1, and < w' < 1. Furthermore, 
the growth condition p.3p implies that C{9) > P as soon as \9\ > R{P) + 3. Indeed, 
in this case \u2 — ui| > 16*1 — 3 > R{P) and right-hand sides of formulas (|5.3p are 
bounded by P from below. Therefore one can set 

e inf{6': C{e) = mmC{d')} > -oo, 9 = sup{6l: C{e) = minC(e'')} < oo, (5.5) 

0' 6' 

where min is attained because C is continuous. 

The set djcgmingi C^p^^^p^^{9') lies on the segment [6,6]. (Indeed, if e.g. 9 < 
6, then C[p„,p^]{9) > C{9) > mine> C{9') > mine' C[p„,p^]{e'), so 9 cannot be- 
long to argmiug/ C[Fq_Fi](^'); ''^ similar conclusion holds ii 6 > 6.) It follows that 
C{p^,P,] (6 - 0) < < C[V^,,^^] (6-1-0). 

°By convexity ^^](^ + 0) < (qF„,Fi](^) - qF„,Fi] (e))/(0 - 6) for all 9 > 6. 
The right-hand side of the latter inequality can be estimated from above by 

L^inf ^(^)-g(Q). (5.6) 

e>e 9-e 



The ratio in the right-hand side takes finite values, so L is finite. This establishes the 
inequality C'^ 
in particular 



inequality C'^^ (6-1-0) < L. The rest of (|5.4p is given by a symmetrical argument; 



e<e Q- 9 

□ 

Definition 5.3. A locally optimal transport plan jg^ is called globally optimal 
if 9o e arg mine C[pg^p^] (9) . 

We can now reduce minimization of (j2.1[) on the unit circle to minimization 
of (|2.2I) on R, which involves the cost function c rather than c: 

Theorem 5.4. The canonical projection tt: M — T establishes a bijection be- 
tween globally optimal transport plans on M x R and transport plans on T x T that 
minimize (|2.1I) . 

Proof. A transport plan 7 on T x T minimizes (|2.ip if it is a projection of a 
transport plan on K x R that locally minimizes the transport cost defined by the cost 
function c{x,y) = minfegzc(x,?/ -I- k) (see introduction; min here is attained because 
of the integer shift invariance and growth conditions (13. 2p . (|3.3I) '). 

Denote S — {(x, y) : c(x, y) = c{x, y)} and observe that the support of the globally 
optimal plan 7^^ lies within S: indeed, if it did not, there would exist a (nonlocal but 
periodic) modification of jg^ bringing some of the mass of each period to S and 
thus reducing the average cost. Therefore jg^ is locally optimal with respect to the 
cost c{x, y) and its projection to T x T minimizes (12.11) . 

Conversely, a minimizing transport plan on T x T can be lifted to R x R in such 
a way that its support lies inside S (translations of arbitrary pieces of support by 
integer increments along x and y axes are allowed because they leave c(a;, y) invariant). 
Therefore its average cost per period cannot be less than that of a globally optimal 
transport plan on R x R. □ 
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6. Fast global transport optimization. In a typical application, such as 
the image processing problem described in the introduction, measures a-nd fii 
come in the form of histograms, i.e., discrete distributions supported on subsets 
X — {xi, 12, ... , Xng} and Y = {yi, y2, . . . , Vm} of the unit circle. These two sets 
may coincide. In what follows we replace X and Y with their lifts to the universal 
cover X and Y and assume that the points of the latter pair of sets are sorted and 
numbered in an increasing order; 

• • • < X-i < Xq = Xna - 1 < < Xi < • • • < a;„„ < 1 < Xno + l = Xi + 1 < . . . , 

• • • < 2/-1 < yo = Vm - l<0<?/i<---< ym < 1 < Vm+i = yi + K ■ . ■ ■ 

Denote masses of these points by ^^{{xi}) = rnf \ = m^p] these are assumed 

to be arbitrary positive real numbers satisfying X]i<i<no — X^KjXm '^f^ ~ ^■ 

6.1. Computation of the average cost and its derivative. Define i{9) as 
the index of min{yj : Ffijjj) > 0} and denote yf = 2/2 = yjie)+i, Vni = 

yj(e)+ni~i- AH the values 

F0(^l),F0(a:2),...,i^0(^no),^^f(yi),^^l'(y2),---,^f(<) (6.1) 

belong to the segment (0, 1]. We now sort these values into an increasing sequence, 
denote its elements by < V(2) < • ■ • < ^'(no+ni) a'^'i ''^(0) 0- Note that for 
each V such that ^(fc-i) < v < V(^k) with 1 < k < no + ni the values = Fq^{v) and 
y(k) = (^i')^^('^) are uniquely defined and belong to X, Y. It is now easy to write an 
expression for the function C^pg^p-^y. 

C[Fo,F^]i9) ^ c(a;(fc),y(fc))(w(fe) (6.2) 

l<A,'<no+ni 

Observe that, as the parameter 9 increases by A9, those 'y(fc) that correspond 
to values Ff decrease by the same increment. Let Ff{yjg) be such a value. As it 
appears in (|6.2p twice, first as V(^k) and then as — in the next term of the sum, it 
will make two contributions to the derivative C'^p^ Fi](^)'- ""^(-^cT^l-^f (j/jo))' J/jo) and 
c{F,~\F^iy,„)),y,„+,) (see fig. Ml top). 

Moreover, there are exceptional values of 6 for which two of the values in (|6.ip 
coincide and their ordering in the sequence (w(fc)) changes. For such values of 6 the 
derivative Cj'^^ p_^^ has different right and left limits, as illustrated in fig. 16.11 bottom: 

C'Fo,F,]iO {c{F^\Ft{yj)),y,+i) - ciF,\F?{y,)),y,)), (6.3) 

l<j<ni 

ClF„.F,]iO + 0) -5](c(Fo-i(i^f(y,) -0),2/,+i) -c(F„-i(Ff(y,) -0),y,)). (6.4) 

If 6 is not exceptional, the value of C'^^^ Fi](^) given by the first of these formulas. 

The function C^p^^p-^] is therefore piecewise affine (see in particular fig. l5.1[ where 
this function is plotted for atomic marginals /io, /ii). Moreover, from the Monge 
condition (|3.ip it follows that Cj^^ Fi](^ — 0) < C'^^^ Fi](^ + ^) exceptional points, 
giving an alternative proof of convexity of C[pg p^^{9) in the discrete case. 

Lemma 6.1. Values of C and its left and right derivatives can be computed for 
any using at most 0{n^ + ui) comparisons and evaluations of c{x,y). 
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A0 



O 



Fn 



Fo 



A0 



F' 



-A9 



o 



F^\Ff{y,))^X V, v,+, 

(a) 9 not exceptional 



Ff 



-AO 



Fo-'iFfivj)) % y,+i 

(b) 9 exceptional, case C'^p^ Fi](^ ~ 




Fo-\F^{y,)~0) y, y,+, 

(c) 9 exceptional, case C'^p^ Fi](^ ~'~ 



Fig. 6.1. Derivation of expressions 116. 3| I. 1 16. 4| I for C'^p^^ ^^j(9±0). Thick lines show fragments 
of complete graphs of Fq, Ff corresponding to jth terms in 116.311 . 116.41 1. thin dashed line (bottom) 
marks the common value of Fq and Ff . Note that X = F~^ {Ff {pj)) is equivalently expressed as 
mf{x: Fo(x)> Ff{y,)}. 



Proof. Sorting the uq + ni values (|6.1I) into an increasing sequence requires ng + 
rii — 1 comparisons (one starts with comparing Fq{xi) and F^{yf) to determine ^{i), 
and after this each of the remaining values is considered once until there remains only 
one value, which is assigned to V(^no+ni) with no further comparison). At the same 
time, pointers to and should be stored. After this preliminary stage, to find 
the values for Cjj?^ j?^] and its one-sided derivatives it suffices to evaluate each of the 
no + ni terms in (|6.2p and to take into account the corresponding contribution of plus 
or minus c{x(^k) , y(k)) to the value of Cj^,^ (0), paying attention to whether the value 
of 9 is exceptional or not. All this can again be done in 0(no + ni) operations. □ 

6.2. Transport optimization algorithm. Fix e > and set L = max{L, L}. 
Recall that L, L, as well as the parameters 9, Q that are used in the algorihtm below, 
are defined by explicit formulas in Lemma [521 and do not depend on measures /io, Mi- 
The minimum of C[pg p-^-^{9) can be found to accuracy e using the following binary 
search technique: 



1. Initially set 0_:— Q and 6 :— Q, where Q, Q are defined in Lemma [5.21 

2. Set 9 ■= lii+d). 

3. Compute qV„,F,](^ " 0), qV„,F,](^ + «)■ 

4. If C'^p^ — 0) < < C'^p^ (0 + 0), then 9 is the required minimum; stop. 
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5. li 6 — d_< e/L, then compute C^p^ p-^] {d_), C^p^^p-^] (9), solve the linear equation 

C[p,.F,]{e) + c[p^^p^]ie + o){e-e) = qp„p,]ie) + c[p^^p^]ie-o){e-0) (6.5) 

for 6, and stop. 

6. Otherwise set := 6* if C[p^,p^] {9 + 0) < 0, or ^ 6* if C[p^,p^] {9 ~ 0) > 0. 

7. Go to step [2 

It follows from inequalities (j5.4p of Lemma 15.21 that the minimizing value of 9 
belongs to the segment [9 , 8] . Therefore at all steps 

C[p„ ,^,] + 0) < < C[p^^p^] (9 - 0) (6.6) 

and the segment [9, 9] contains the minimum of C. 

Step [5] requires some comments. By convexity, —L < Cj^^ {9 ± 0) < L for 

aU e < < e, i.e., |qV^^^](6l ± 0)1 < L at all steps. When 9 - 9 < e/L, this 

bound ensures that for any 9' in [9,9] the minimal value of C is within e/L ■ L = e 
from C[pg,p-^]{9'). If there is a single exceptional value of 9 in that interval, then it 
is located precisely at the solution of (j6.5l) and must be a minimum of C because 
of (j6.6p , so the final value of 9 is the exact solution; otherwise it is an approximation 
with guaranteed accuracy. 

The final value of 9 will certainly be exact when masses of all atoms are rational 
numbers having the least common denominator M and e < l/M. Indeed, in this case 
any interval [9_, 9] of length e can contain at most one exceptional value of 6. 

Since at each iteration the interval [6_, 9] is halved, step [5] will be achieved in 
0(log2((8 — 6)/(e/L))) iterations. By Lemma l6T] each instance of step[3](and equa- 
tion ()6.5p ) takes O{no + ni) operations. Thus we obtain the following result. 

Theorem 6.2. The above binary search algorithm takes 0{{no + ni) log(l/e)) 
comparisons and evaluations ofc{x,y) to terminate. The final value of 9 is within e/L 
from the global minimum, and C[pg^p^^{9) < ming C[pg p-^-^{9) + e. When all masses 

mf'\ vri'p are rational with the least common denominator M, initializing the algo- 
rithm with e = 1/2M leads to an exact solution in 0((no + ni) log Af) operations. 

6.3. Experiments. We tested experimentally the estimates of Theorem l6.2l for 
time complexity as a function of parameters of the problem. The average computing 
time of the algorithm for different values of uq + ui and e is illustrated in fig. 16.21 
These results have been obtained using the following procedure. For each value of 
no and ni, points {xi } and {yi, ^2, ■ • ■ , yni}, which constitute the sup- 

port of distributions fiQ and fii, are drawn independently from the uniform distribu- 
tion on [0,1] and sorted. The masses rnf'\ i = 1 . . . no and ni^^\ j = l...ni are 
then drawn from the uniform distribution and normalized such that X)i<i<no = 

J2i<j<ni "^^P — 1- Finally the transport cost is minimized for c(x,y) = \x — y\. The 
code used to produce this figure is available online at the web site of the OTARIE 
project http : //www .mccme . ru/~cLnsobol/otarie/ software .html. 

In the first experiment the value of e was set to 10^^", the algorithm was run 
10 times for each pair (no, ni) with 1 < no, ni < 100, and the computing times were 
averaged. In the second experiment (noi^^i) was fixed at (10,10) and the average 
computing time was similarly computed for different values of e. The averaged com- 
puting times for the two experiments are plotted in fig. 16.21 Observe the manifest 
linear dependence of computing time on no -I- ni and log e. 



18 



J. DELON, J. SALOMON AND A. SOBOLEVSKI 




100 120 140 160 180 200 




Average computing time, sec., vs no + ni for (b) Average computing time, sec, vs log^g e for 
: 10-1". no = ni = 10. 



Fig. 6.2. Average computing time of the algorithm for different values of (no +ni) and logjQ e. 
The experiment was performed on a PC with a 3.00 GHz processor. 




(a) Top: F. Maliavin, Whirlwind (1916). 

(b) Right: P. Puvis de Chavanne, Jeunes fllles au 
bord de la mer (1879). 




(c) Puvis' hues transported to Mahavin's canvas. 
Fig. 6.3. Optimal matching of the hue component of color. 



FAST TRANSPORT OPTIMIZATION ON THE CIRCLE 
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The next figure is a concrete, if not entireiy serious, iilustration of optimal matcli- 
ing in the case of distributions on the "color circle." 

Recall that in the HSL (Hue, Saturation, and Lightness) color model, the color 
space is represented in cylindrical coordinates. The polar angle corresponds to the 
hue, or the degree to which a color can be described as similar to or different from 
other colors (as opposed to difference in saturation or lightness between shades of the 
same color). 

We chose two famous paintings, one Russian and one French, whose highly dif- 
ferent coloring is characteristic of the two painters, the Expressionist Filipp Maliavin 
(1869-1940) and the Symbolist Pierre Puvis de Chavannes (1824-1898). An optimal 
matching of the hue distributions according to the linear cost c{x,y) = \x ~ was 
used to substitute hues of the first painting with the corresponding hues of the second 
one while preserving the original values of saturation and brightness. In spite of the 
drastic change in coloring, the optimality of matching ensures that warm and cold 
colors retain their quality and the overall change of aspect does not feel arbitrary or 
artificial. 

7. Related algorithmic work. Fast algorithms for the transportation problem 
on the circle, with the Euclidean distance |a; — yj as a cost, have been proposed in 
a number of works. Karp and Li |13| consider an unbalanced matching, where the 
total mass of the two histograms are not equal and elements of the smaller mass 
have to be optimally matched to a subset of elements of the larger mass. A balanced 
optimal matching problem has later been considered independently by Werman et al 
[20] : clearly, the balanced problem can always be treated as a particular case of the 
unbalanced one. In both of these works O(nlogn) algorithms are obtained for the 
case where all points have unit mass. 

Aggarwal et al [T] present an algorithm improving Karp and Li's results for an 
unbalanced transportation problem on the circle with general integer weights and the 
same cost function \x — y\. They also consider a general cost function c{x,y) that 
satisfies the Monge condition and an additional condition of bitonicity: for each x, 
the function c{x,y) is nonincreasing in y for y < yo{x) and nondecreasing in y for 
y > yo{x). Note that this rules out the circular case. The second algorithm of [T] 
is designed for bitonic Monge costs and runs in 0(n log M) time for an unbalanced 
transportation problem with integer weights on the line, where M is the total weight 
of the matched mass and n is the number of points in the larger histogram. 

The algorithm proposed in the present article only applies to the balanced problem 
for a Monge cost. However it does not involve bitonicity and is therefore applicable 
on the circle, where it achieves the same 0(n log M) time as the second algorithm 
of [1] if all weights are integer multiples of 1/M. Although our theory is developed 
for the case of costs satisfying a strict inequality in the Monge condition, it can be 
checked that the discrete algorithm works for the case c{x, y) = \x — y\, which can be 
treated as a limit of \x — y\^, A > 1, as A ^ 1 17J. 

Finally we note that results of [T] were extended in a different direction by Mc- 
Cann [16], who provides, again in the balanced setting, a generalization of their first 
algorithm to the case of a general cost of the concave type on the open line. This 
case is opposite to Monge costs and requires completely different tools. Indeed, for a 
strictly concave cost such as c{x, y) = \J\x — y\ the notion of locally optimal trans- 
port plan on the universal cover does not make sense: concave costs favor long-haul 
transport over local rearrangements, destroying local finiteness. 
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